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12,  13  &  14  June  2007 


Thomas  A.  Donnelly,  Ph.D.  (ECBC), 
Erin  E.  Shelly  (ECBC),  and 
Daniel  P.  Cinotti  (SAIC) 


UNCLASSIFIED/UNLIMITED 


wily  4Js^  Design^fiElcperiments  (DOE) 
Methods  with  Simulaiion  Experiments? 


Quicker  answers,  lower  costs, 
solve  bigger  problems 

O  Obtain  a  fast  and  cheap  surrogate  “meta-model”  of 
the  simulation 

■  can  more  rapidly  answer  “what  if?”  questions 

■  do  sensitivity  analysis 

O  By  running  efficient  subsets  of  aii  possibie 

combinations,  one  can  -  for  the  same  resources  and 
constraints  -  soive  bigger  probiems 

O  Be  as  cost  effective  as  possibie  and  run  no  more 
triais  than  are  needed  to  get  a  usefui  answer 
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•  Demonstrated  how  Design  of  Experiments  (DOE) 
can  be  used  to  sequentialiy  run  groups  of 
simulation  trials  to  obtain  better  and  better 
meta-models  of  the  simulation  model 

•  When  control  variables  are  all  continuous  and 
response  variable  is  NON-stochastic,  then 
“Smoothing”  designs  can  be  used  to  efficiently 
produce  a  meta-model  of  a  simulation  that  is 
made  up  of  a  complex  series  of  physical  models 


UNCLASSIFIED/UNLIMITED 
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ZJ  ,  UNQLASSIFIED/UNLIMITED 

- - iTvY<^-Types-6f  Designs  for 

T4/v6  Types  or  Meta^iiMioc^eling  of  Simulations 


•  “Traditional”  designs  for  polynomial  modeling 
with  categorical  and  continuous  variables 

O  Designs  can  be  sequentially  constructed  to  support 
increasingly  complex  models 

O  Featured  example  reanalyzes  a  simulation  case  matrix  in  which 
all  648  combinations  of  variable  settings  were  originally  run 

•  “Smoothing”  designs  for  use  with  continuous 
variables  AND  non-stochastic  responses 

O  Though  little  used,  these  designs  are  a  more  efficient  alternative 
to  traditional  designs  and  exploit  “Kriging”  regression  analysis 


UNCLASSIFIED/UNLIMITED 
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•  If  a  “textbook”  fractional-factorial,  orthogonal 
array  or  response-surface  design  is  available, 
then  use  it. 


•  Textbooks  and  web  site  catalogs  do  not  always 
contain  designs  for  categorical  variables  with: 

O  all  combinations  of  mixed  numbers  of  ieveis  (e.g.  3,  4,  5,  and  21) 

O  iarge  numbers  of  levels  for  variables  (e.g.  5+) 

•  Algebraic  (Orthogonal  Array)  and  algorithmic  (D-optimal) 
computer  generated  designs  can  often  be  used 

O  Orthogonal  Arrays  are  good  at  yielding  analysis  with  “clean” 
(unconfounded)  estimates  of  the  “main  effects” 

O  D-optimal  designs  are  good  for  adding  on  the  fewest  additionai  trials  to 
support  higher  order  “interaction”  terms  in  the  modei 

UNCLASSIFIED/UNLIMITED 
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is6  Matri)^  (TBM  Bulk)  & 

^LIMLTED 

Example  Dosage  Plot 

. I  " 

i 

I 

las 

Used  ih  $tudy  of  the 

dbserved  Response 

^‘ProbabilRyoTCasualty”  (PCAS) 

Variable 

#  Levels 

Levels 

Agent  Codes'' 

6 

A,  N,  T,  H,  R,  Y  (categorical) 

Season 

3 

Winter,  Summer,  Spring/Fall  (categorical) 

Time  of  Attack 

3 

0500,  1200,  2200  Local  Time  (continuous) 

No.  of  TBMs  & 
Spread  Radius^ 

2 

1  TBM  &  1  m,  2  TBMs  &  1000  m  (categorical) 

Mass^>^  (relative) 

3 

1.00,  1.57,  2.00  (continuous) 

Height  of  Bursts 

2 

0,  10  m  (continuous) 

Total  Cases 

648 

1. 

2. 

3. 

4. 

5. 


Dropped  “Q”  ■  it  had  smallest  effect  &  6  levels 
allowed  for  use  of  a  smaller  Orthogonal  Array 
Spread  Radius  paired  with  No.  of  TBMs 
Mass  (with  3  levels)  replaced  Source  Strength 
(with  2  levels) 

Mass  is  nested  in  Agent 

Data  was  available  for  Height  of  10  m 


Target  Area 


See  Presentation  by  Dan  Cinotti 
in  WG-2,  NBC  Defense, 
1100  on  Thursday 
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•  Because  a  different  set  of  mass  values  were 
used  for  each  agent,  the  variable  Mass  is 
“nested”  within  the  variable  Agent 

•  The  response  Probability  of  Casualty  (PCAS), 
which  is  bounded  within  the  range  (0,  1),  was 
transformed  using  2*Arcsin((PCAS)''^^)  which 
maps  the  range  (0,  1 )  to  the  range  (-<»,  +<») 

O  This  made  the  error  fit  the  usual  regression  assumption  of 
being  normally  distributed 

O  This  also  prevented  our  regression  from  predicting  values  and 
limits  that  were  above  1.0  and  physically  impossible 
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■ 

1 

1 

Foul 

r  Stage  Design  Sequence 

stage  1 


Stage  2  Stage  3  Stage  4 


36  Total 
Simulations 


108  Total 
Simulations 


324  Total  ALL  648 

Simulations  Simulations 


Design  1, 36  trials 


Main  effects  only 
for  ALL  variables 


Stage  1  effects 
plus  all  2-way 
interactions 


Design  1, 36  trials 
Design  2, 72  trials 


Design  3, 216  trials 

Stage  2  effects 
plus  all  3-way 
interactions 


5.6%  of  648 


16.7%  of  648 


50%  of  648 


324  trials  in  Design  4  used  as  checkpoints  for  Designs  1,  2  &  3 

UNCLASSIFIED/UNLIMITED 


Design  1, 36  trials 
Design  2, 72  trials 


Design  3, 216  trials 

Stage  3  effects 
plus  ALL 
remaining  4-way, 
5-way  and  6-way 
interactions 

Design  4, 324  trials 

NOTE:  Length  of  this 
green  box  should  be 
longer  than  shown 
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1 

i 

T^b 

-Pre 

dictJohS'  oiF  PCAS 

vs 

1 

i>r 

Plot  (Continuous) 
2nd  Order  Model 

Interaction  -  PCAS 


0.80^ 

0.75 


1 

1.5714 

2 

MunCnLS  = 

c 

) 

A 

o.9r? 

0.959 

0.974 

H 

0.916 

0.935 

0.944 

N 

0.966 

0.986 

0.992 

R 

0.833 

0.868 

0.384 

T 

0.730 

0.772 

0.793 

Y 

0.930 

0.986 

0.989 

Interaction  -  PCAS 


1.0  1.2  1.4  1.6  1.8  2.0 

Mass 


1 

1.5714 

2 

MunCnt S  = 

Q 

) 

A 

0.977 

0.991 

0.995 

H 

0.947 

0.960 

0.966 

N 

0.993 

0.998 

0.999 

R 

0.395 

0.920 

0.931 

T 

0.302 

0.837 

0.354 

Y 

0.990 

0.993 

0.995 

Mass 


Season  =  F 
Time  =  12.0 
HOR  =  0.0 

MunCnt Spi 


2 


A 

N 

T 
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77 I  ,  n ! I 

Predlcti'^s 

r  UNCLASSIFipAUNLIM^ED 

(w/95%  Pred.  Limits)  of  PCAS  vs. 

Nested  Mas 

»s  arpd  MdhCdtlSpread  for  1-way, 

reduced  2-way  and  reduced  3-way  models 

1-way  Model,  Highlighted 
Prediction  is  0.802  ±  0.030 
Based  on  fitting  ^  trials 


2-way  Model,  Highlighted 
Prediction  is  0.803  ±  0.008 
Based  on  fitting  108  trials 
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3-way  Model,  Highlighted 
Prediction  is  0.802  ±  0.000 
Based  on  fitting  324  trials 
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:  II  I  -  I  UNCLASSIFIED/UNLIMITED 

Percent  Off  Target  for  3|24  PCAS  Checkpoint 
Pr^d^timis  With  jl -Way,  |2-Way  and  3-Way  Modeis 

“How  Gofodns  Good  Enough?” 


1-way  Model 
Fit  to  36  Trials  in 
Stage  1  Design 


Common  Scale  range  for 
plots  is  from  -4%  to  4% 


PCTOFF1WAY 


Worst  Case  =  3.7% 
Half  of  Cases  <  0.37% 


Reduced  2-way  Model 
Fit  to  36  +  72  Trials  in 
Stage  2  Design 


Worst  Case  =  -0.93% 
Half  of  Cases  <0.11% 


Reduced  3-way  Model 
Fit  to  36  +  72  +  216  Trials  in 
Stage  3  Design 


Worst  Case  =  0.008% 
Half  of  Cases  <  0.001  % 


UNCLASSIFIED/UNLIMITED 
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^  -  unglassih^/unlimTted 

Seminal  Bapen4>n  |"Smoothing*”  DOE 
for  Computer  Experiments 


•  Sacks,  J.,  Welch,  W.J.,  Mitchell,  T.J.  and  Wynn, 
H.P.  (1989).  “Design  and  Analysis  of  Computer 
Experiments.^^  Statistical  Science  4.  409-423 

O  First  textbook  appeared  in  2003  and  has  the  same  name 

O  A  good  source  for  up-to-date  information  is  the  Simulation 

Experiments  &  Efficient  Designs  (SEED)  Center  for  Data  Farming  at 
http://harvest.nps.edu 

*Smoothing  is  an  alternate  name  sometimes  used  for  designs  for  computer 
experiments  because  it  is  a  good  description  of  the  end  result  of  the  analysis. 
Another  name  that  sometimes  appears  is  “space-filling”  designs  because 
trials  are  spread  somewhat  uniformly  throughout  the  test  volume. 
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•  From  the  traditional  experimental  design  point  of 
view  the  Smoothing  designs  -  for  the  same 
number  of  trials  -  do  not  enclose  as  large  a 
volume  of  the  design  space.  This  is  intentional. 

•  Rather  than  emphasizing  high  leverage  trials 
(“corners”)  for  a  simple  polynomial  model,  these 
designs  “spread”  their  trials  more  uniformly 
through  the  space  to  better  capture  the  local 
complexities  of  the  simulation  model. 

•  Analysis  employs  “Kriging”  method  originally 
developed  for  geo-spatial  regression 


UNCLASSIFIED/UNLIMITED 
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Prl>cess 


UNGLASSIFJIEDZUMLIMITED 


|itiniizauoni4>f-Modelecl  Industrial 


u4ii^gXonliputer  Experiments 


•  Data  is  generated  by  a  simulation  consisting  of 
a  series  of  physical/chemical  models  each 
feeding  its  result  into  the  next. 

•  Industrial  examples  include: 

O  Chemical  plant 
OAircraft  engines 
ODeep  ocean  oil  production 
O Semiconductor  fabrication  line 
OAluminum  can  extruder 

Ran  51  “designed”  simulation  trials,  analyzed  data,  determined 
optimal  factor  settings,  checked  optimum  with  a  simulation  trial 
(they  agreed),  built  1  real  machine  for  $500,000  and  made  real  cans 
-  the  performance  was  “dead  on” 

•  DoD  examples  include  M&S  like  the  ECBC 
Chem-Bio  Sim  Suite,  SOES  Smoke  Model,  etc. 


UNCLASSIFIED/UNLIMITED 
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C.3  Examples 


The  following  examples  (iemonstrate  many  ])()ssil)le  uses  {)f  PErK.  The  re¬ 
sponses  for  these  examples  are  based  on  the  Brmiiii  fiuiction.  The  Braniti 
function  is  the  real-valued  function  of  tw'o  variables 


FIGURE  C.l.  The  Branin  function  on  [—5, 10]  x  [0, 15] 
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Example  Latin  Hypercube  Design  and 
Data  Caiculated  with  Branin  Function 


1 

7-75 

6 

35.80951 

2 

1 

3.75 

14.86287 

3 

10 

8.25 

31.41880 

4 

4.75 

4.5 

19.87899 

5 

2.5 

15 

141.88566 

6 

-3.5 

2.25 

99.43335 

7 

3.25 

0 

3.88973 

8 

■5 

6.75 

97.47380 

9 

-4.25 

12.75 

6.27060 

10 

6.25 

1.5 

19.85914 

11 

8.5 

11.25 

95.50587 

12 

7 

14.25 

181.74214 

13 

-0.5 

0.75 

49.39445 

14 

-2 

5.25 

23.13762 

15 

0.25 

10.5 

43.09524 

16 

9.25 

3 

2.82392 

17 

-2.75 

9.75 

3.61474 

18 

5.5 

9 

75.79100 

19 

4 

12 

104.11175 

20 

-1.25 

13.5 

43.33586 

21 

1.75 

7.5 

23.39797 
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C.3  Examples 


T\w  following  Dxaiijplos  (lcnK)tist  rafc  many  [jossiblc^  nsos  of  PErK.  The  re¬ 
sponses  for  these  (examples  are  on  Ihe  Ilnniht  JiniriioiL  The  Braniii 

function  is  llie  real- valued  function  of  t  wo  varial>les 


Trial  X1  X2  Y 


1 

7.75 

6 

35.80951 

2 

1 

3.75 

14.86287 

3 

10 

8.25 

31.41880 

4 

4.75 

4.5 

19.87899 

5 

2.5 

15 

141.88566 

6 

■3.5 

2.25 

99.43335 

7 

3.25 

0 

3.88973 

8 

■5 

6.75 

97.47380 

9 

-4.25 

12.75 

6.27060 

10 

6.25 

1.5 

19.85914 

11 

8.5 

11.25 

95.50587 

12 

7 

14.25 

181.74214 

13 

-0.5 

0.75 

49.39445 

14 

■2 

5.25 

23.13762 

15 

0.25 

10.5 

43.09524 

16 

9.25 

3 

2.82392 

17 

-2.75 

9.75 

3.61474 

18 

5.5 

9 

75.79100 

19 

4 

12 

104.11175 

20 

-1.25 

13.5 

43.33586 

21 

1.75 

7.5 

23.39797 

UNCLASS 


100  150  200  250  300  350 
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FIGURE  C.l.  The  Branin  function  on  [-5, 10]  x  [0, 15] 


4(X) 


2(X) 


Plot  from  software  of  “Kriging” 


0  -6 


fit 


^0 
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Comp^ring|  Suitjaces  for  Increasingly  Complex 
Po|yiioi|nia1s|nt~to  Da^a  ^rpm  the  Branin  Function 


60() 


4(M) 


2(K) 


Linear 

y  =  ao  +  +  82X2 


600 


40t) 


200 


Interaction 
y  =  Linear  +  a  ^2^  1^2 


600 


4(X) 


7^ 


0  A 


0  A 


600 


4(X) 


2(K) 


Quadratic 
y  =  Interaction 

+  +  322X2^ 


)  Partial  Cubic 

6(X 

Cubic 

600 

y  =  Quadratic 

y  =  Partial  Cubic 

)  +  a^^2^i^^2  ^122^1^2^ 

4(K 

)  +  a^y^Xy^  +  a 222^2^ 

40() 

2(X 

200 

ric'^ 

“Smoothing” 
or  “Kriging”  Fit 


0  Jb 


The  full  cubic  model  appears  to  closely  approximate  the  Branin  function, 
but  still  cannot  represent  the  ripples  seen  in  the  fit  using  Kriging  method. 
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UNCLASSIFIEO/UNLIMITED 


CB  Suite  Sii^oot hi ng  DOE 
Exampie  with  1 0  Variabies 


•  Branin  function  example  is  trivial.  With  2  control 
variables  the  full  cubic  model  has  10  terms. 

•  The  following  example  has  10  control  variables. 
(Full  cubic  model  has  166  terms!) 

•  Three  different  Smoothing  designs  are  used: 

1. 17-trial  Latin  Hypercube  (LHC)  design 

2.  33-trial  Nearly  Orthogonal  Latin  Hypercube  (NOLH) 
design  (see  SEED  web  site  at  http://harvest.nps.edu) 

3.  50-trial  Orthogonal  Array  (OA)  design. 

•  Smoothing  design  trials  combine  in  such  a  way 
as  to  fall  into  5  of  6  Pasquill  Atmospheric 
Stability  regions  within  the  VLSTRACK  model 
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Uill 
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Atmospheric  Stability  Classes  & 
li/le^eorologi^al  Conditions  That  Define  Them 


stability  Class  Definition 

A 

very  unstable 

B 

unstable 

C 

slightly  unstable 

D 

neutral 

E 

slightly  stable 

F 

stable 

Key  point  is  that  VLSTRACK 
models  each  class  a  bit  differently 
and  we  want  to  create  a  single 
meta-model  of  all  classes  together 


Surface  Wind  Speed 

Daytime  Incoming  Soiar  Radiation 

Nighttime  Cloud  Cover 

m/s 

mi/hr 

strong 

Moderate 

Slight 

>  50% 

<  50% 

<  2 

<  5 

A 

A-  B 

B 

E 

F 

2  to  3 

5  to  7 

A-  B 

B 

C 

E 

F 

3  to  5 

7  to  11 

B 

B  -  C 

C 

D 

E 

5  to  6 

11  to  13 

C 

C  -  D 

D 

D 

D 

>  6 

>  13 

C 

D 

D 

D 

D 

Note:  Class  D  applies  to  heavily  overcast  skies, 

at  any  windspeed  day  or  night 

TABLES  SOURCE  :  http://en.wikipedia.Org/wiki/Air_pollution_dispersion_terminology#_note-7#_note-7 


ORIGINAL  SOURCE:  Pasquill,  F.  (1961).  The  estimation  of  the  dispersion  ofwindborne  material, 
The  Meteorological  Magazine,  vol  90,  No.  1063,  pp  33-49. 
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IIKVMViVUIdldiTilJ.^llIJilUdil 


CB  Simulation  Suite  Architecture 


Threat  Delivery 


Hazard  Environment 


Real-time  Sensors 
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CB  Sim  Suite  is  a  set  of  distributed  simulation  tools  designed  to  represent  all  aspects  of  CB 
passive  defense  on  the  tactical  battle  field  for  application  to  analysis,  testing,  and  training. 
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-irliJcioi  s-in  ]2-D-fi>r 
iLl^mlooithin^”  besigns  o 


/UNLIMITED 


Different  10-Variable 
biz^  17,  33  &  50  Trials 


17-trial 

Latin  HyperCube  (LHC) 


33-trial 
Nearly  Orthogonal 
Latin  Hypercube  (NOLH) 


50-trial 

Orthogonal  Array  (OA) 


UNCLASSIFIED/UNLIMITED 
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_  UNCLASSIFIED/UNLIMITED 

50-trial  Orthogonal  Array 
with  5  Levels  per  Variable 


Showing  first  17  of  50  trials  in  one  “space-filling”  design 
out  of  =  9,765,625  possible  combinations  of  variable  settings 


UNCLASSIFIED/UNLIMITED 
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LIMITED 


krlging  lAiialysis  of  a  oingle  Simulation  - 
Concentratibn  vs.  Latitude,  Longitude  &  Time 


Cloud  release  point  is 
1 0  km  west  of  1 0  km 
X  10  km  grid  of  72 
identical  entities 


Wind  speed  is  5.3  m/s 
Wind  direction  is  278° 
from  north 


Concentration<S> 


25 


25 


Lon(coded) 

25  25 


25 


25 


28  min.  32  min.  36  min.  40  min.  44  min. 

UNCLASSIFIED/UNLIMITED 


2C 


5 


48  min. 

25 


Wind  direction 


kriging 


¥- 


Ahalysis  of  1 7\ 


Simulations 
atiohs 


U^ii)g  [17  0bs 
Max  Dosage  vs.  8  Variables 


Dosage<S> 


1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0 


time  wrt  sunset  =  360. i 
Temp  =  21.0 
Humidity  =  50.0 
CloudCover  =  0.620 
Amount  Agent  =  200.0 
Log(cluration)  =  1.0 


Wind_speed 

Wind_spe=3.80  Wind_dir=270.00 

Value 

311.82 


1200 


Time  wrt  Sunset  =  360 
Temperature  =  21 
Humidity  =  50 
Cloud  Cover  =  0.62 
Amount_Agent  =  200 
Logio(Duration)  =  1 


\N\nd_speed 


UNCLASSIFIED/UNLIMITED 
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rigmg 


UNCLASSIFIED/UNLIMITED 


Ai{ialysis  of  ilTi  L^C  Simulations 
□sing  120^  Observation^  =  17  X  72  -  15 
Max  Dosage  v^.  10  Variables 


Dosage<S> 


time  wrt  sunset  =  360. 
Temp  =  21 .0 

Humidity  =  50.0 
Cloud  Cover  =  0.620 
Amount  Agent  =  200.0 
Log(duration)  =  1.0 
Latitude(coded)  =  17.0 
Longitude(coded)  =  17.C 


Wind_speed 

Wind_spe=3.80  Wind_dir=270.00 

Value 

404.61 _ 

1200  1200  1200 


Cloud  Cover  =  0.26  0.38  0.50 


1200 


Dosage<S> 


80  3 


Time  wrt  Sunset  =  360 
Temperature  =  21 
Humidity  =  50 
Cloud  Cover  =  0.62 
Amount_Agent  =  200 
Logio(Duration)  =  1 
Latitude  =  17 
Longitude  =  17 


UNCLASSIFIED/UNLIMITED 
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In  FutiJre  Will  jSliow  %  Off  T^rgpt  for  200  Checkpoint 
Predi^tioris  ^ith  Varidu^  Sknoothing  Designs 
“Suspect  -  Never  Get  Something  for  Nothing” 


Kriging  Model  Kriging  Model  Kriging  Model 

Fit  to  17  Trials  in  Fit  to  17 +33  Trials  in  Fit  to  17  +  33  +  50  Trials  in 

LHC  LHC  +  NOLH  LHC  +  NOLH  +  OA 


Common  Scale  range  for 
plots  is  from  -4%  to  4% 


PCTOFF1WAY 


Worst  Case  =  ?% 
Half  of  Cases  <  ?% 


-4  -2^0  2 

PCTOFF2WAY 

Worst  Case  =  ?% 
Half  of  Cases  <  ?% 


PCTOFF3WAY 


Worst  Case  =  ?% 
Half  of  Cases  <  ?% 
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10-Variable  Meta-Model  Predicting  Concentration 


100 


Off-Axis  Variable  Settings 

Time  wrt  Sunset  =  360 
Wind  Speed  =  3.8 
Wind  Direction  =  270 
Humidity  =  50 
Cloud  Cover  =  0.50 
Logio(Duration)  =1.0 
Latitude  (coded)  =17 
Longitude  (coded)  =17 

NOTE:  This  is  a  plot  of  Kriging  regression  of 
the  100  integers  between  0  and  99  randomly 
assigned  to  100  smoothing  design  trials. 

The  “noise”  has  been  fit  perfectly! 

This  is  why  one  should  only  use  this 
technique  with  non-random  data! 
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•  Demonstrated  how  Design  of  Experiments 
(DOE)  can  be  used  to  sequentiaiiy  run  groups 
of  simuiation  triais  to  obtain  better  and  better 
meta-modeis  of  the  simuiation  modei 

•  When  controi  variabies  are  aii  continuous  and 
response  variabie  is  NON-stochastic,  then 
“Smoothing”  designs  can  be  used  to 
efficientiy  produce  a  meta-modei  of  a 
simulation  that  is  made  up  of  a  complex 
series  of  physical  models 
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